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Abstract 

The transport properties of an ionic model for liquid silica p| at high tem- 
peratures and pressure are investigated using molecular dynamics simulations. 
With increasing pressure, a clear change from "strong" to "fragile" behaviour 
(according to Angell's classification of glass-forming liquids) is observed, al- 
beit only on the small viscosity range that can be explored in MD simulations. 
This change is related to structural changes, from an almost perfect four- fold 
coordination to an imperfect five or six-fold coordination. 
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More than 10 years ago, Angell [§] proposed a classification of glass forming liquids ac- 
cording to their "fragility", a concept that has since then proved extremely useful for our 
understanding of these systems. Broadly speaking, " strong" glass formers are characterized 
by weak discontinuities of their thermodynamic properties at the glass transition, and an 
almost Arrhenius behaviour of their transport properties (viscosity) as a function of tem- 
perature. The corresponding activation energy E a can be assigned to a typical microscopic 
"event' that controls transport properties, such as breaking of a "bond". The archetype of 
"strong liquids" is S1O2, with a perfectly Arrhenius behaviour of the viscosity in the range 
1500 < T < 2500. The activation energy is in that case about 70000K, and is supposed to 
correspond to the energy for breaking an Si — O bond. 

Fragile liquids, on the other hand, display a strongly non Arrhenius behaviour of their 
transport properties. The flow in these systems can not be traced back to a single microscopic 
event, but is of a much more collective nature. Typical fragile liquids are organic or ionic 
systems with a much weaker local organization than in "strong" liquids, it is now widely 
accepted that at least for relatively low viscosities, the transport properties in these systems 
are well described by the so called "mode-coupling" theories ||. 

It has long been suspected, and also shown in some pioneering MD simulations [|J, 
that the transport properties of S1O2 could exhibit an unusual behaviour as a function of 
pressure, in particular, a nonmonotonous, "waterlike" behaviour of the diffusion coefficient 
is expected due to the tetrahedral network structure of molten Si0 2 ■ It is also expected ]5],|2| 
that the disruption of this network by an applied pressure will cause a change in behaviour 
from "strong" to "fragile". 

In this work, we present a systematic Molecular Dynamics (MD) study of the transport 
properties in a purely ionic model of Si02 ||1| , which is an improvment of an earlier model 
by Tsuneyuki and coworkers ||. Both models have now been studied in much detail for 
their static properties 0J|,|H],|n],[| . Transport properties at zero pressure have also been 



briefly investigated in [13], and their pressure dependence was considered in [I5|. 

The system we consider is made up of 216 Si02 units (in some cases this number was 



increased to 512 in order to check the size dependance of the results). The equations of 
motion are integrated using the Verlet algorithm with temperature control [12| . Four series of 
runs at four different densities (2.2, 3.1, 3.3 and 4.2 g/cc) were carried out. For each density, 
the sample is equilibrated at high temperature (5000K), then cooled stepwise at an average 
cooling rate of about 10 12 K.s -1 . At each new temperature, the system is equilibrated during 
0.25 to 0.5 nanosecond. As usual f|, this equilibration time limits the temperature range for 
which the system can be considered to be at thermodynamic equilibrium. By analogy with 
the situation in real systems, the lowest temperature at which the system is equilibrated for 
an equilibration time of 0.5 nanosecond defines the "computer glass transition temperature" 
for our simulation. In terms of viscosities, this corresponds to viscosities close to 1. Pa.s.. 

The self diffusion coefficients of Si and O, D$i and Do, were computed from the mean 
squared displacement of the corresponding ions. The shear viscosity 77 was computed from 
the stress-stress autocorrelation function using the usual Kubo formula ||. The results for 
r) and 1/Dgi are shown in figures 1 and 2, in the form of an Arrhenius plot. As usual, the 
error bars on rj (estimated from the contributions of the 3 components of the stress tensor) 
are large, typically 100% for the larger viscosities. The results, however, are unambiguous. 
At the lowest density (2.2 g/cc), the Arrhenius plot is almost a perfect straight line, with an 
activation energy about 40000 Kelvins. As previously noted by Delia Valle and Andersen 
0, this activation energy is smaller than the experimentally observed one. Most likely, this 
constitutes an inadequacy of the purely ionic model. Another interesting possibility, however, 
is that the difference is due to the different temperature range investigated. Experiments 
are done at low temperatures, while simulation investigates high temperatures. In fact, a 
slight curvature of the Arrhenius plot is perceptible if the last point, corresponding to a 
temperature of 3200K, is included. This would mean that the activation "energy" changes 
with temperature. Unfortunately, the error bar on this point is large, due to the very small 
diffusivity. High temperature experiments on Si02 would be useful to clarify this issue. As 
the density is increased, a decrease in rj and 1/D$i is observed, while the Arrhenius plot 
displays a noticeable curvature. At the higher density, the viscosity has increased close to 
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its low density value, but the Arrhenius plot is clearly bent, in a way characteristic of fragile 
liquids. The nonmonotonous variation of D$i with pressure is illustrated in figure 3 for a 
fixed temperature (T = 4000-fT). As pressure is increased, Dgi increases by almost an order 
of magnitude, with a maximum in the 10-15GPa range. This effect was already observed 
in |T3j (with a different model) and [|T5(] . The size effect reported in the latter paper is 
not observable in our results, as seen from the two points in figure 3 obtained with a 1536 
particles system. 

Diffusion and vicosity are often related by introducing a "Stokes Einstein diameter" 
dsE of the atoms, defined as dsE = hsT /ZirqD. This quantity, shown (for the Si ions) in 
figure 4 for the four densities investigated, is seen to be fairly independent of temperature. 
The scatter observed at low temperature for the lower density is probably due to the large 
error bars on the corresponding viscosities. At the lowest density, d S E is large (3 to 4 A), 
but becomes much smaller, and almost independent of density(l to 1.5 A), for the higher 
densities. This illustrates the very peculiar character of low density S1O2, with a perfectly 
tetrahedral coordination. The evolution of the local structure with pressure is illustrated in 
figure 5, which shows the proportion of four, five and six fold coordinated Si as a function of 
pressure at T = 4000K. In the standard picture of diffusion in Si0 2 , five-fold coordinated 
Si are often described as "defects" of the tetrahedral network that favour diffusion. Figure 
5 shows that the proportion of these "defects" increases very rapidly with pressure, and 
that in the vicinity of the diffusion maximum a large majority of Si ions are indeed five-fold 
coordinated. Hence it seems unlikely that a " defect based" picture of Si0 2 can be useful in 
describing the pressure dependance of its properties, except perhaps at very low pressures. 
Figure 5 also shows the evolution of coordination with pressure, from fourfold at P = 
to fivefold in the 10-15 GPa range and sixfold above 20 GPa. The low pressure four fold 
coordination, however, is very special in the sense that it is alost perfect, with more than 
99% of the ions. At higher pressures, the coordination is much less well defined. 

Finally, the influence of the local structure on self diffusion is illustrated in figure 6. 
The mean squared displacement of Si ions with different initial coordinations (3, 4, or 5) 
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is plotted as a function of time. For long times, these quantities become parallel straight 
lines. However at short times, the difference between the different curves is a measure of 
the influence of coordination on the diffusion. This role is very clear in low density S1O2, 
where three or five fold coordinated Si (defects) diffuse initially much faster than the four- 
coordinated Si. At a density of 2.8g/cc, (corresponding to about 6Gpa in figure 3 and 5), 
the role of initial coordination has become almost negligeable. The overshoot in the mean 
squared displacement that exists at low densities and disappears at higher densities was 
noticed in |L£| and associated to the existence of a "boson peak". This overshoot, however, 
might be associated to a finite size effect |HJ . 

In summary, our results indicate that the disruption of the tetrahedral network of low 
density S1O2 by the application of pressure very rapidly induces a change from strong to 
fragile behaviour. It must be remembered, however, that the change observed here corre- 
sponds to a viscosity variation of less than three decades. Experimentally, the same kind 
of change could be expected to tkae place in a much larger viscosity domain, covering more 
than ten decades. In the vicinity of the diffusivity maximum, the behaviour is already that of 
a fragile liquid. The role of five or three fold coordinated "defects" in the diffusion becomes 
negligeable before the diffusivity maximum is reached. 

These findings corroborate those obtained with other ionic models of S2O2 ||16|| . They 
confirm the richness of these simple models. Beyond the fact that they provide a relatively 
accurate description of silica, a great interest of these models lies in the fact that they involve 
simple pair potentials, for which accurate theoretical methods have been devised. They could 
therefore provide a useful benchmark for theoretical studies of strong/fragile transitions. 
Finally, the connection between these results and experiment could be made in the high 
temperature-high pressure range (typically T ~ 2500K, P ~ 20GPa. This range can be 
explored experimentally, and corresponds (according to the model) to diffusion constants 
and viscosities measurable on the simulation time scale. A comparison with experimental 
data would allow a better evaluation of the usefulness of such models in the prediction of 
dynamical properties. 
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FIGURE CAPTIONS 
Figure 1: Arrhenius plot of l/D Si at four different densities. Open dots: p = 2.2g/cc. 
Open stars: p = 3.3g/cc Stars: p = 3.7g/cc. Squares: p = 4.2g/cc. The lines are a guide to 
the eye. 

Figure 2: Same as figure 1, for the viscosity. The symbols are as in figure 1. 

Figure 3: Pressure dependance of the silicon diffuion constant at T = 4000 A. The 
points corresponds to densities p = 2.2, 2.5, 2.8, 3.1, 3.4, 3.7, 3.9, 4. 2g/cc. For p = 2.5g/cc 
and p = 3Ag/cc two different system sizes (N=256 and N=512 Si0 2 units ) were considered. 

Figure 4: Stokes-Einstein diameter of the silicon ions as a function of temperature, for 
the same four densities, as in figures 1 and 2. Symbols as in figure 1 and 2. 

Figure 5: Fraction of 4, 5, and 6 coordinate silicon ions as a function of pressure, for 
T = 4000 A. The densities are the same as in figure 3.. Inset: the equation of state of the 
simulated system at T=4000K 

Figure 6: Mean squared displacement of silicon ions that are initially 3,4, and 5 fold 
coordinated. Dots: fourfold coordinance. Long dashes: fivefold coordinance. Short dashes: 
threefold coordinance. From top to bottom: p = 2.2g/cc, p = 2.8g/cc, p = 3.1g/cc. T = 
4000A. 
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